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ABSTRACT 

The thermal evolution of neutron stars (NSs) is investigated by coupling with 
the evolution of r-mode instability that is described by a second order model. The 
heating effect due to shear viscous damping of the r-modes enables us to understand 
the high temperature of two young pulsars (i.e., PSR B0531+21 and RX J0822- 
4300) in the framework of the simple npe NS model, without superfluidity or exotic 
particles. Moreover, the light curves predicted by the model within an acceptable 
parameter regime may probably cover all of the young and middle-aged pulsars in 
the Ig — \gt panel, and an artificially strong p superfluidity invoked in some early 
works is not needed here. Additionally, by considering the radiative viscous damping 
of the r-modes, a surprising extra coohng effect is found, which can even exceed the 
heating effect sometimes although plays an ignorable role in the thermal history. 
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1 INTRODUCTION 

The composition of neutron star (NS) interior is still poorly known due to the uncertainties of 
nuclear physics. It has been hoped that comparing theoretical thermal evolution of NSs with ther- 
mal emission data from observations could yield information about their internal properties. Our 
knowledge of the cooling history of a NS has been improving as we were refining the physical 
ingredients that play a key role on the thermal evolution of NSs. Page et al. (2004) proposed the 
minimal cooling model and suggested that four basic physical inputs should be considered in the 
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simulation of NS themial evolution, namely, the equation of state, superfluid properties of relevant 
components, envelope composition and stellar mass. Aguilera et al. (2008) argued that the mag- 
netic field is another basic physical ingredient should not be ignored. In fact, heating mechanisms 
should not be ignored, too. 

During the evolution of NSs, heating mechanisms may be present and dramatically change the 
thermal evolution of the star. Several heating mechanisms have been extensively discussed, for 
example, rotochemical heating (Reisenegger 1995; Fernandez & Reisenegger 2005), vortex creep 
heating (Umeda et al. 1995), Joule heating (Aguilera et al. 2008), heating due to the hardron-quark 
transition (Kang & Zheng 2007) and heating due to r-mode damping (Zheng et al. 2006; Yu et al. 
2009). 

i?-modes in a perfect fluid star with arbitrary rotation arise due to the action of the Coriolis 
force with positive feedback (Andersson 1998; Friedman & Morsink 1998), succumbing to grav- 
itational radiation-driven Chandrasekhar-Friedman-Schutz instability. However, in a realistic star, 
the r-mode evolution is determined by the competition between damping effect due to viscous 
dissipation and the destabilizing effect due to gravitational radiation. Based on the conservation 
of angular momentum, a phenomenological model describing the r-mode evolution was proposed 
by Owen et al. (1998) and improved by Ho & Lai (2000). In this original version of the model, 
an unbounded growth could lead the modes to an unphysical regime because nonlinear effects 
are ignored. As an important nonlinear effect, differential rotation induced by r-modes was first 
studied by RezzoUa et al. (2000, 2001) and confirmed by some numerical studies (Stergioulas & 
Font 2001; Lindblom et al. 2001). Sa (2004) solved the fluid equations within nonlinear theory up 
to the second order in the mode amplitude and described the differential rotation analytically. By 
extending the r-mode evolution model of Owen et al. (1998) to this nonlinear case, Sa & Tome 
(2005, 2006) obtained a saturation amplitude of r-modes self-consistently. Yu et al. (2009) studied 
the long-term spin and thermal evolution of isolated NSs under the influence of the differential ro- 
tation. They found that the stars can keep nearly a constant temperature for over a thousand years 
since the differential rotation can significantly prolong the duration of r-modes. 

In this paper, we study the thermal evolution of isolated NSs with differential rotation induced 
by r-modes using a realistic equation of state (EOS). The same EOS is used by Kaminker et 
al. (2001) in the simulation of NS cooling, but their interpretation of the observations requires a 
strong p superfluidity. As pointed by Tsuruta et al. (2002) and Blaschke et al. (2004), the strong p 
superfluidity contradict the microphysical calculations in case of high enough proton concentration 
to permit nucleon direct Urea process. The inclusion of r-mode dissipation would probably resolve 
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this contradiction. In this work, we also take into account the recently realized radiative viscosity 
(Sa'd et al. 2009). 

After the introduction in this section. Section 2 presents the formulism of r-mode evolution and 
the spin evolution of NSs up to second order in the mode amplitude. Section 3 displays the equation 
of thermal evolution and the heating (or cooling) term due to viscous dissipation of r-modes. The 
physical inputs and our results are given in Section 4, and the conclusions and discussions are 
presented in the last section. 



2 THE i?-MODE EVOLUTION AND THE SPIN EVOLUTION OF NSS TO 
SECOND-ORDER 

For a rotating barotropic Newtonian star, the r-mode solutions of perturbed fluid equations can be 
found in spherical coordinates (r, 6, <p) at first order in a as (Lindblom et al. 1998), 

6^^V = 0, (1) 
J^V = aaCil{^^^ \m'-^esm(l(P + a)t), (2) 

5^1^^* = QrQC//(^j' \in'-^6lcosecos(/0 + a>O, (3) 
and at second order in a as (Sa 2004) 

^(2)/ = 5^%^ = 0, (4) 

WOAr^-' sin^-i 6, (5) 

where a represents the amplitude of the oscillation, R and Q. are the radius and angular velocity of 
the unperturbed star, oj = -Q.(l+2)(l- !)/(/+ 1), Q = (21- 1)! ! V(2/ + l)/[2;r(2/)!/(/ + 1)], A and 
are two constants determined by the initial condition. For simplicity, Sa & Tome (2005) suggested 
N = 21- I and redefined A by introducing a new free parameter ^ as A = ^KCjf(l + 1)R^~^K For 
the most unstable / = 2 r-mode of primary interest to us, the second-order solution shows 
a differential rotation of the star induced by the r-mode oscillation, i.e., large scale drifts of fluid 
elements along stellar latitudes. 

Using 6^^^v' and 6^^^v', the corresponding Lagrangian displacements and can be derived 
and then the physical angular momentum of the / = 2 r-mode can be calculated up to the second 
order in a as (Sa 2004; Sa & Tome 2005) 

J, = + /2) = ^^E±^a^jMR^a^ (6) 
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where / = 1.635 x 10"^ and 

J''^ = -Jpd,e^d,^fUv'V,^f')dV, (7) 

+v'(Vi^f + V,^f)]dV. (8) 

Meanwhile, following Owen et al. (1998) and Sa (2004), we further express the energy of the 1 = 2 
r-mode by 

= _ = ^^El^a^jMR^a\ (9) 

3 2 

When K = -2, 7^^^ vanishes and the expressions of Jr and Er return to their canonical forms (Owen 
et al. 1998), in other words, the differential rotation disappears. Both the physical angular momen- 
tum and energy of r-modes are increased by gravitational radiation back reaction and decreased 
by viscous damping, which yields 
dJj_ _ 2J^_2Jr^ 

dt Tg Ty 

dEy ^ IE,- lEf (]-\\ 

dt Tg Tv 

where Tg is the growth timescale due to gravitational- wave emission, t„ = {t~1 + t^I + t'^IY^ is 
the damping timescale due to viscous dissipation; Tv,„Ti,v and t„ are the timescales of the shear 
viscous damping, bulk viscous damping and radiative viscous damping (the radiative viscous will 
interpreted latter), respectively. 

From equation (flOl) and (fTTI) . it can be seen that the r-modes are unstable if {t~^ -t'^Y^ > 0. In 
this case, a small perturbation would lead to a non-ignorable growth of the modes. The competition 
between the gravitational destabilizing effect that is dependent on Q. and the T-dependent viscous 
damping effect determines an instability window in the T - Q. plane. 

For a normal NS with a strong magnetic field (~ 10^°~^^ G), besides the braking effect due to 
gravitational radiation, the spindown of the star resulting from magnetic dipole radiation should 
also be taken into account. So, we ought to write the decrease of the total angular momentum of 
the star as (Owen et al. 1998; Ho & Lai 2000; Sa & Tome 2005) 
dJ 3a^JMR^Q IQ 



(12) 



where = 1.35 x 10^5j^^(Q/ y/nGpY^s, is the magnetic braking timescale and / = IMR^ with 
/ = 0.261 is the moment of inertial of the star. Due to the r-mode oscillation, the total angular 
momentum of the star could be separated into two parts, i.e., J = IQ. + Jr. Then, Eqs. (fTOl) and (fT2)) 
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yield 

da 
dt 



l + -{K + 2)Qa' 



a 



1 + \^iAK + 5)Qa' 



a ^ a 



dQ. 8 ,Q 2 ,Q Q 
— = --(K + 2)Qa^- + -{4K + 5)Qa' , 



(13) 
(14) 



where Q = 3 J/21 = 0.094. 



3 THERMAL EVOLUTION OF NSS 

The equation of thermal evolution can be written as (Yakovlev et al. 1999; Yakovlev & Pethick 
2004) 
dT 

Cv— = -Ly -Ly+ H,. - ALy, (15) 
dt 

where Cy is the total stellar heat capacity. In NSs composed of simple npe matter, the electrons 
constitute an almost ideal, strongly degenerate, ultra-relativistic gas; neutrons and protons consti- 
tute a non-relativistic strong non-ideal Fermi liquid. 

Meanwhile, Ly is the luminosity of neutrinos generated in numerous reactions in the interiors 
of neutron stars. The main processes we used are nucleon direct Urea, nucleon modified Urea and 
nucleon bremsstrahlung. If the proton and electron Fermi momenta are too small compared with 
neutron Fermi momenta, the nucleon direct Urea process is forbidden because it is impossible 
to satisfy conservation of momentum (Lattimer et al. 1991). Under typical conditions, one finds 
that the ratio of the number density of protons to that of nucleons must exceed about 0.11 for the 
process to be allowed. 

Ly is the surface photon luminosity given by 

Ly = 4nR^o-T^, (16) 

here cr is the Stefan-Boltzmann constant and is the surface temperature.The relation between 
and the internal NS temperature T is taken from Potekhin et al. (1997), supposed the outer heat 
blanketing NS envelope is made of ion and neglecting the effects of surface magnetic fields. We 
also note that the effective surface temperature detected by a distant observer is T'^ = ^Jl - Rg/R, 
where Rg is the gravitational stellar radius. 

(Hy - ALy) is the energy per unit time induced by viscous dissipation of r-modes. As we can 
see from equation ([TT]) . one part of the oscillation energy of r-modes is converted into heat energy 
through shear viscous damping and bulk viscous damping(//v = 2Er{T'~^l + t^I)), and the other 
part is converted into neutrino emissivity through radiative viscous damping(ALv = 2£',-t~!). As a 
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result 



1 1 1 

H, - AL, = 2Er{— + ), 



(17) 



^.vv "^bv "^rv 

using equation (|23l) below, we get 



1 1 1 
H,-ALy = 2Eri -— ). 



(18) 



Hence, if Tgy < iT^y, we have (//,, - ALy) > 0, the star would be heated by viscous dissipation of 
r-modes ; and if Tw > 2tm,, we have {H^ - ALy) < 0, viscous dissipation causes an extra cooling 
of the star. 

4 PHYSICAL INPUTS AND RESULTS 

4.1 NS EOS 

In the simulation of thermal evolution, we employ the simplest possible nuclear composition, 
namely neutrons (n), protons {p) and electrons (e). We adopt a moderately stiff equation of state 
(EOS) of this matter proposed by Prakash et al. (1988) (their model I with the compression modu- 
lus of saturated nuclear matter K = 240MeV). The maximum mass of this model is M = 1.911 Mq, 
and the direct Urea process is forbidden at M < Mq = 1.358Mo. 

4.2 The r-mode timescales 

The calculation of r-mode timescales are very complicated, and until now there are no calculations 
based on realistic EOS. The following timescales (for 1 = 2 r-modes) we employed in this paper 
are obtained with a polytropic equation of state as ;? = kp^ for NSs, with k chosen so that the 
mass and radius of the star are M = lAM^ and R = 12.53 km. In the following equations, the 
conventions = T/IO^ and d = Q./ ^JnGp are used. 

The gravitational radiation timescale is (Owen et al. 1998) 



The shear viscous damping timescale (due to the neutron-neutron scattering) is (Owen et al. 
1998) 



In the case of npe matter where only the modified Urea process is relevant, the bulk viscous 
damping timescale is (Andersson & Kokkotas 2001) 



Tg = 3.26Qr^s. 



(19) 



T,, = 2.52 X lO^Th. 



(20) 
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Thv(MUrca) = 1.20 x lO^^T~^Cl-^s, 



(21) 



and in the case of direct Urea proeess (Owen et al. 1998) 



rUDUrca) = 6.99 x lO%^Cl-^s. 



(22) 



Recently, Sa'd et al. (2009) first demonstrated that there exists a new mechanism for damping 
the energy of stellar oscillations, namely the radiative viscous dissipation. Urea processes con- 
tribute to the damping of density perturbations not only by converting energy into heat via bulk 
viscosity, but also by converting it into an increase of the neutrino emissivity via radiative vis- 
cosity. They found the radiative viscosity coefficient is 1.5 times larger than the bulk viscosity 
coefficient. Thus, the damping time scale of radiative viscosity is 

Trv = ^'^bv (23) 

Since Tbv(MUrca) is about three orders larger than TtviDUrca), it is natural that for the two 
kinds of NSs the direct Urea process is permitted or not, the r-mode instability windows (a window 
determined by r^^ - r^T^ = in the T - Q. plane) are quite different. In contrast, the time evolution 
behavior of a and Q. for these two different kinds of NSs are similar. We can see this from eqs. (fT3l) 
and (fT4l) (the following argument is based on the relation Tg « Tm, which is correct unless the 
magnetic fields of NS is extremely large). For a nascent NS with T ~ 10'°K and ~ f y/^Gp, 
we can easily find Tg « t^v and Tg « Thv This means the terms of Tg in eqs. (fT3]) and (fT4l) 
are decisive in the early evolution of NS for both cases r^, = Ti„(DUrca) and r^, = Tbv(MUrca). 
More exactly, fig. 3 in Yu et al. (2009) (they adopted Tfo,, = TtviDUrca)) indicates that not only for 
a nascent NS but also for NS of phase I, II and III, Tg dominates. Meanwhile, during the following 
phase IV and V in their figure, t^v « Thv, Tg and t^v control the evolution of NSs. Thus, we 
conclude that ti„ plays an ignorable role during the evolution of NSs. 

4.3 The results 

We calculated Equations ([13] ). (fT?l) and (fTSl ) numerically, taking the initial temperature Tq = 10'°K, 
the initial r-mode amplitude ao = 10"^, the initial angular velocity fio = | y/nGp and the magnetic 
field B = 10^2 Q 

Fig.l shows the evolution of a, Q. and T of a IAMq NS and K = 1000. We can see from Fig.l 
(c) that the NS core can keep high temperature for 4.6 x 10^(or 10^^^)years. As illustrated by Yu 
et al. (2009), during the early part of the r-mode evolution, the rotation energy of the star ijIO?') 
is converted into the oscillation energy, the internal energy, and the energy of gravitational waves. 
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Nevertheless, during the late part, the energy deposited in the r-modes would be released gradually 
via heating the star and accelerating the stellar rotation. 

Fig.2 (the parameters are the same as Fig.l) displays the evolution curves of (H.^, - ALy) due to 
viscous dissipation and (-Ly - Ly). The evolution of {Hy - ALy) can be divided into three stages: 
(I) Igt < -4.61. In this phase, the relation r^v > 2Tby is fulfilled, and the viscous dissipation of r- 
modes results in an extra cooling to the thermal evolution of NS. This extra cooling never plays an 
important role during the evolution of the NS because it presents in the early part of NS evolution 
where the star interior is still very hot, and it is too small comparing with the neutrino luminosity 
Ly. Moreover, it's easy to understand that the same results can be reached to NS of small mass 
(M < I.358M0) where the direct Urea process is forbidden. (2) -4.61 < Igt < 3.65. In this phase, 
Tsv < 2Tby, the r-mode energy is dissipated mainly by shear viscosity and the star is heated. (3) 
Igt > 3.65. In this phase, (Hy - ALy) = 0, since the r-modes disappear in this stage (see Fig. 1(a)). 

Fig. 3 shows the surface temperature T'^ evolution of M = 1.3Mq, M = 1.365Mo and M = 
IAMq NS with different K. Note that the direct Urea process is forbidden at M < Md = 1.358Mo 
in our EOS model. We can see that all curves taking into account r-mode dissipation can explain 
two young and hot pulsar data (PSR B053 1-^21 and RX J0822-4300) when taking proper value of 
K(K = 100 for I.3M0 and M = 1.365M0 , K = 1000 for I.4M0 ). Maybe the I.3M0 (K = 100) 
curve looks too high for the explanation of the two data. While, it should be emphasized we 
don't expect our model to explain the observation data precisely; because this model still have 
some uncertainties, such as the realistic non-linear evolution of r-mode instability is not clear 
and the timescales we used are not based on the specified EOS. Nonetheless, it is certain that the 
dissipation of the oscillations can provide large enough heat to raise the temperature of young 
neutron stars, and this isn't a negligible effect. 

5. CONCLUSIONS AND DISCUSSIONS 

We have studied the thermal evolution of NSs, considering the influence of r-mode instability to 
second order. For the first time we take into account the radiative viscosity, and find the radiative 
viscous dissipation of r-modes results in an extra cooling of NSs; while in former studies it is 
taken for granted that viscous damping would lead to the heating of stars. However, we find that 

this extra cooling can be well neglected in the thermal evolution history of NSs. 

On the other hand, the NS is heated due to shear viscous damping of r-modes, and it can keep 
a high temperature for several thousand years, even tens of thousands of years. This enables us 
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to explain two young and hot pulsar data (PSR B0531+21 and RX J0822-4300) with NS model 
composed of only npe matter, without superfluidity or exotic particles. In contrast, under the same 
NS EOS, Kaminker et al. (2001) explained these data by the inclusion of strong p superfluidity (the 
maximum critical temperature T^p'^Sx lO^K). The superfluidity they employed is too strong and 
it was doubted by many works (Tsuruta et al. 2002; Blaschke et al. 2004); because Takatsuka & 
Tamagaki (1997) showed, through careful microphysical calculations, that for neutron matter with 
such high proton concentration as to permit nucleon direct Urea process, the superfluid critical 
temperature should be extremely low, about several xlO'K. However, if consider a wider value 
range of NS mass and K, our light curves may probably cover all of the young and middle-aged 
thermal emission data, and the artificially strong p superfluidity invoked in Kaminker et al. (2001) 
is no longer needed . Therefore, the explanation of observation data maybe doesn't contradict 
the calculation of microphysics. Of course, we don't expect our result to fit the observation data 
accurately, but future studies using an improved model, or even including other heating effect 
would further improve our results. 
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Figure 1. Evolution curves of a, n, and T of a 1.4Mo NS with magnetic field B = lO'^ G and = 1000. 
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Figure 2. Evolution curves of (Hy - ALy) due to viscous dissipation (tiiick solid line) and (-Ly - Lj) (dashed line) with the same parameters as 
Fig.l. The thin soUd one is the zero luminosity line. 
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Figure 3. Observational data on the surface temperatures of NSs (Yakovlev et al. 2008) compared with theoretical coohng curves. The dot, dashed 
and solid curves correspond to M = I.SMq, M = 1.365Mo and M = IAMq, respectively. The thick hues are calculated without the r-mode 
dissipation effect, and the thin lines refer to K = 10, 100, 1000 form left to right. For all curves, B = lO'^G. 



